Effective non-additive pair potential for lock-and-key interacting particles: the role of 

the limited valence. 
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Theoretical studies of self-assembly processes and condensed phases in colloidal systems are often 
based on effective inter-particle potentials. Here we show that developing an effective potential for 
particles interacting with a limited number of "lock-and-key" selective bonds (due to the specificity 
of bio- molecular interactions) requires — beside the non-sphericity of the potential — a (many body) 
constraint that prevent multiple bonding on the same site. We show the importance of retaining both 
valence and bond-selectivity by developing, as a case study, a simple effective potential describing 
the interaction between colloidal particles coated by four single-strand DNA chains. 

PACS numbers: 61.20Gy, 61.20.Ja, 82.70.Dd 



Description of multi-scale phenomena in self-assembling 
systems and understanding of their collective behavior re- 
quires a coarse graining process in which the irrelevant 
degrees of freedom are integrated out in favor of an ef- 
fective re- normalized interaction potential [1]. Colloidal 
systems are often described by this type of approach, in 
which the solvent and its chemical-physical properties are 
encoded into parameters of the resulting effective inter- 
colloid particle potential [2j]. Thermodynamic (and often 
also dynamic) bulk behavior of several complex particles, 
including polymers, star-polymers, dendrimers, micelles, 
telecheric polymers 0, 0, HI, & 0, Hj] has been interpreted 
on the basis of effective pair-wise potentials, function only 
of the distance between particle centers. 

Very recently, design of new materials has started to 
capitalize on the possibilities offered by complementary 
"lock-and-key" selective interactions [9], aiming to achieve 
a much better control on the resulting bonding pattern and 
on the mechanical properties of the self-assembled bulk 
material. In this Letter we show that effective potentials 
for these new particles must account for both the valence 
(maximum number of bonds) and for the selectivity of the 
bonding process. While at the microscopic level speci- 
ficity in the bonding interactions is pair-wise additive — 
it builds on steric incompatibilities [10], e.g. on excluded 
volume effects which prevent the possibility of multiple 
interactions at the same bonding site — we find that the 
conservation of this important feature at the level of the 
coarse-grained effective potential may require the break- 
down of the pair-wise additive approximation. 

As a working example, we discuss the development of 
an effective potential for a macromolecule composed by a 
core particle decorated by a small number of short sin- 
gle strands (ss) DNA [11] to show (i) the importance of 
properly accounting for valence and bond selectivity; (ii) 
the need for a non-pair-wise additive effective potential. 
The full model is graphically explained in Fig. [T] (and 
the chosen length and energy units), and more in detail 
Ref. [l2|, [l3[. A semi-rigid chain of eight particles, model- 
ing a 8-base ss-DNA, is attached to each of the four parti- 
cles composing a central tetrahedral core. Each of the 32 



base particles carries an additional interaction site (labeled 
A,T,G or C) which can bind only with complementary site 
types (only A-T or G-C binding are allowed), resulting in 
a total of 68 interaction sites per macromolecule. The cho- 
sen base sequence is palindromic, so that all eight bases 
can be simultaneously paired when ss-arms of different 
macromolecules approach in the correct orientation. All 
site-site interactions are pair-wise additive. The phase di- 
agram of this model shows a region of gas-liquid phase 
separation, but only at small number densities p. The 
liquid region is characterized by the presence of a four- 
coordinated network, whose dynamics progressively slows 
down on cooling, parallel to the formation of an open gel- 
like structure [131 ]. 

To evaluate the effective potential we simulate, using 
Monte Carlo methods, a system composed of two macro- 
molecules in a box for several temperatures T, covering 
the T-region where ss-DNA pairs. From a large ensem- 
ble of equilibrium configurations we extract the properties 
of the degrees of freedom to be retained in the coarse- 
grained model: (i) the normalized probability h r of find- 
ing the two particles at relative center-to-center distance 
r (i.e. in such a way that -p- J 47rh r r 2 dr = 1) (ii) the 
normalized probabilities h$ 1 and h$ 2 of the angle 0\ and 
O2 between f and the closest ss-DNA arm on each parti- 
cle (so that - JJ 7 hesinOdO = 1). The variables r, 0\ and 
62 are defined in Fig. [TJ The T-dependence of h r (r) and 
he(0i) is shown in Fig.[2j On cooling, angular and spatial 
correlations build up, reflecting the formation of a stable 
double-strand configuration, with a characteristic center- 
to-center distance of the order of the DNA arm length 
(8.5 < r < 10.5) and a preferential orientation for linear 
bonding (#1 « 0). At the highest and lowest investigated 
T, h r and he describe, to a good approximation, the fully 
non-bonded (nb) and fully bonded (b) state, respectively. 
The non-bonded high T state is characterized by essen- 
tially unstructured hg b and h™ h except for the r-region 
close to the origin where excluded volume interactions play 
a role. The function h b e is well described by a gaussian 

distribution h h e = 2 e ~^ (where the normalization 
constant is A = 22.32 and the best-fit a 2 = 0.052 rad 2 ). 
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FIG. 1: Drawing of two macromolecules in a bonded con- 
figuration. Each macromolecule is composed by a tetrahedral 
hub (four tetrahedrally-oriented orange spheres) to which four 
single strands of eight bases each (grey spheres) are anchored. 
Each base is decorated by a short-range attractive site (small 
spheres) encoding the base type (A,T,G,C). Base complemen- 
tary is enforced by selective attractive interactions between A 
and T and between G and C sites. In the process of devel- 
oping a coarse-grained model, we associate a rigid tetrahedron 
(oriented according to the central core of the particle) to each 
macromolecule and focus on the center-to-center distance r and 
on the two angles 0\ and #2, defined in the figure. In this Letter 
we use as unit of length the diameter of the larger spheres (or- 
ange or grey) and as a unit of energy the A-T (or the equivalent 
G-C) binding energy. Number density is defined as the total 
number of large spheres (36 for each macromolecule) divided 
by the volume. 



A similar angular distribution has been recently chosen to 
model self-assembly of patchy particle [14J. The distribu- 
tion h b is essentially localized only between 8.5 < r < 10.5, 
corresponding to the bonded distance. We have not been 
able to provide a simple analytic best-fit expression for h b 
and h™ b and, in the following, we have used their numer- 
ical values. The sharp separation in configuration space 
between the bonded and the non-bonded states (due to 
the localized nature of the interaction and to its speci- 
ficity) suggests to describe the distributions at interme- 
diate T as a linear combination of the b and nb distri- 
butions, weighted respectively by the probability of being 
bonded (p&(T)) or non-bonded (1 — pb{T)). Lines in Fig. [2] 
show that both h r and he can indeed be described as lin- 
ear combinations (with the same Pb(T) coefficient) of the 
b and nb distributions. Moreover, as shown in the in- 
set of Fig. [21 the resulting T-dependence of Pb(T) is very 
well fitted by the two-state expected theoretical expression 

_ AU-TAS 

Pb(T) = 1/(1 + e k B T ), where AU and AS measure 
the change in energy and entropy associated to the forma- 
tion of a double strand. The important result is that the 
T-dependence decouples from the r and dependence and 
it becomes possible to write the probability of finding the 
two macromolecules at relative distance r and angles 0\ 
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FIG. 2: Results from the simulation of two macromolecules in 
a box of side L = 25: center-to-center probability h r (r) (a) and 
angle 0\ distribution he(0i) {he{02) = he(0i) for symmetry). 
Symbols are numerical results for different T. Lines are the 
corresponding expressions provided by Eq. [2] after integration 
over 6\ and O2 in panel (a) or integration over r and O2 in panel 
(b). The inset shows the T dependence of pb and its fit to a 
two-state model (see text) with AE = 4.94 and AS/ks = 52 
(kB is the Boltzmann constant). 



and 02 as 



P{r,0 1 ,0 2 ,T) = 
p h (T) P b (r, e u 2 ) + [1 - PbiT)) P nb (r, U 2 ). 



(1) 



We approximate P b {r,0 1 ,0 2 ) = ftj(r)hj(8i)ftj(fi 2 ) and 
analogously P nb for the non-bonded case. In this approx- 
imation, correlation between distances and angles within 
the b (or the nb) state are neglected, since both P b and 
P nb are written as a product of their arguments. We have 
tested that this is indeed a good approximation. Correla- 
tion between angles and distances arising from the differ- 
ences between b and nb states (small 0\ and 2 are found 
when r is within the bonding distance) are retained in P, 
which indeed can not be factorized. 



Thus the four functions h b , h b e , h™ b 



hf and p h (T) fully 

specify P(r,0i,0 2 ,T) and make it possible to model the 
macromolecule as a particle decorated by four tetrahe- 
drally oriented arms and evaluate an effective potential 
between these particles. The simplest approach is to ne- 
glect the angular correlation and develop a spherical (s) 
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FIG. 3: (a) Shape of the non- spherical effective potential 
/3Vns(r) for given angles 0± = and 62 = 0, and (b) the spher- 
ically averaged effective potential /3V s (r), for different temper- 
atures. 



effective potential V s (r, T) , integrating out the angular de- 
grees of freedom as 

8V(rT) In Pb(?)h b r (r) + [1 - Pb (T)]h?(r) 
PVs[r,l)- m pb{T)h t {roo) + pb{T)]h n b{roo) ■ W 

Here is any distance larger than the range of interac- 
tion of the two molecules (so that V s (r, T) = at large 
distances). A more accurate approach requires retaining 
information on the angular orientations, defining a non- 
spherical (ns) angular dependent potential V ns as 



, P(r,O u 2 ,T) 
11 P(rooA,W ] 



By construction, both V s and V ns accurately describe the 
radial distribution function (and, in the case of V nsi also 
the angular distribution function) of the system in the 
limit of vanishing density at all T. The radial shape of 
the two resulting potentials is shown in Fig. El V ns has a 
deeper minimum than V 3 , to compensate for the reduction 
of the solid angle associated to bonding. 

To assess the ability of the two potentials in reproducing 
the bulk behavior of the original model and to estimate the 
importance of the directional interactions, we compare the 
phase diagram of the two potentials, with the (previously 
studied) behavior of the orig inal model [II Eg. Both V s 
and V ns perform very badly, predicting a very wide region 
of liquid-gas instability, which includes the region where 
the tetrahedral network of bonded particles is known to be 
the stable phase (see Fig. H]). The unstable gas-liquid re- 
gion is even wider in the case of V ns , due to its deeper min- 
imum. The disagreement is so large that no state points 
can be chosen to compare the original and the effective 
models in the region where a network develops, since V s 
and V ns always generate a phase-separated configuration. 
The failure in reproducing the bulk behavior can be as- 



FIG. 4: Phase diagram of the model for colloidal parti- 
cles coated by four single-strand DNA chains. The dashed 
(red) area approximately indicates the region where phase- 
separation is expected for the model studied in Ref. [13[ (es- 
timates based on the extrapolated behavior of the low wave- 
vector structure factor). Filled squares, outside the dashed 
area, denote the stable state points studied in Ref. [13]. The 
graph also shows the boundary between homogeneous and 
phase-separated states for V s and V ns . Note that most of the 
filled squares are located in a region of the phase-diagram where 
phase separation is observed in V s or V ns . The filled circles in- 
dicate the gas-liquid coexistence for the valency constrained 
effective model. 



cribed to the absence, in both V s and V ns , of the lock-and- 
key character of the interaction. In the real model, once 
a ss-DNA is linked to a complementary strand of a differ- 
ent particles, it is not available for further bonding. The 
specificity of the interaction implies that each ss-DNA can 
not be simultaneously bonded to more than one other arm. 
This important rule is not enforced neither by V s (which 
allows for a large number of bonded neighbors) nor by V ns 
since the width of the angular part of V ns does not pre- 
vent multiple bonding at the same arm direction. We are 
forced to abandon the pair- wise additivity approximation 
and develop a new potential which enforces the bond se- 
lectivity and the valence (so that each particle can not be 
engaged in more than four bonds). To this aim, we com- 
plement V ns with an ad-hoc rule imposing that each arm 
can be engaged only in one bond, very much in the same 
spirit as the maximum valence models [15] . From a com- 
putational point of view, during the simulation, we retain 
a list of all pairs of bonded arms and make sure that each 
arm is involved at most in one bonded interaction. The re- 
sulting effective valency-constrained non-spherical model 
is able to generate homogenous equilibrium states in the 
expected region of the phase diagram. This opens the pos- 
sibility to compare the predictions of the effective model 
with those of the original model at the same T and 
p. A comparison in real and in reciprocal space is shown 
in Fig. [5l The limited-valence effective model is able to 
reproduce with satisfactory accuracy the structure of the 
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FIG. 5: Radial distribution functions g(r) for the original 
model (points, from Ref. [13]) and the valency-constrained ef- 
fective potential (lines) at p — 0.03 (a) and p — 0.05 (b). Inset: 
structure factors S(q) at the same state points. Note the in- 
crease in S(q) at small q when p = 0.03, signaling the proximity 
of the limit of stability of the liquid phase. 



tetrahedral network: both the peculiar ratio between the 
position of the first and second peak in the center-center 
radial distribution function g(r) as well as the development 
of a pre-peak in the corresponding structure factor S(q) 
on cooling are well captured by the effective potential. To 
locate precisely the gas-liquid coexistence for the limited- 
valence effective potential we perform umbrella sampling 
grand canonical Monte Carlo simulations [16] (and his- 
togram re- weighting techniques The resulting gas- 
liquid coexistence curve is shown in Fig. [H As compared 
to the unconstrained V s and V nsi the region of gas-liquid 
instability is now confined to small p, opening up a win- 
dow of densities (above the liquid coexistence branch) in 
which equilibrium liquid states can be accessed. 



In summary, we have shown that in the process of devel- 
oping effective potentials for particles with low valence and 
with selective interactions (as in functionalized colloids) it 
is mandatory to account for the specificity of the inter- 
action which intrinsically limit the maximum number of 
bonds that can be formed. The present study shows that 
colloidal particles with specific interactions constitute an- 
other class of colloidal systems for which pair additivity 
approximation breaks down in the effective potential, in 
addition to the recently reported case of charge-stabilized 
colloidal crystal [l8[. The importance of retaining, be- 
side the directionality, the valence of the original model is 
clearly evidenced by the fact that valence strongly controls 
the width of the gas-liquid coexistence [TH, HO] • 
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